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1. TECHNICAL OBJECTIVES 


1. Optimal Utilization of Laser and VLBI Observations for Reference 
Frames for Geodynamics (Grant NSG 5265) 

2. Optimal Utilization of Satellite-Borne Laser Ranging System (Grant NSG 5265) 

3. Geodetic Utilization of NAVSTAR Geodetic Positioning System (Grant NSG 5265) 

4. Utilization of Range Difference Observations in Geodynamics (Contract 
NAS 5-25888) 


2. ACTIVITIES 

2.1 Comparison of Data from Project MERIT Short Campaign 


2.11 Prediction of Polar Motion 

Based on the analysis of polar motion behavior, we found the pos- 
sibility of predicting polar motion for a long time interval (1-2 years in 
advance) with sufficient accuracy. We take the best estimated Chandler 
period as constant, use six years of data to estimate the amplitudes, 
phases and ellipticity (only for annual) of Chandler and annual motion. 

These estimated parameters are then used to predict the next year's (or 
next two years') polar motion. In making the prediction, we also take the 
linear trend into consideration. 

The data used for prediction are BIH, IPMS, and DMA, with the time 
interval from 1968 to 1980 (DMA data from 1972 to 1980). We predicted 
polar motion one yar in advance, compare the predicted polar coordinates 
with the real observed ones (smoothed); the mean rms of the differences 
(predicted minus observed) is about 0'.'02. The differences of the relative 
polar motion are much smaller. For a time interval of 20-30 days, the rms 
of relative polar motion differences is about O'.'Ol (30 cm) through the whole 
year. Compared with the best available VLBI results (from 1977 to 1980), the 
rms of predicted - observed is 0'.'013; and the relative rms (with time inter- 
val less than or equal to two months) is O'.'OOS; here the VLBI observed data 
is unsmoothed. 
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The predicted polar motion can be used for geodetic purposes. It 
also seems that the accuracy of polar motion prediction is high enough for 
any purposes that require the real time polar motion value, including the 
control of spaceships. 

Theoretically, we can say that about 90 percent of polar motion is 
composed of the stable, predictable Chandler and annual terms. 

We also analyzed the error sources in polar motion prediction. In 
doing so, we found that in DMA data there are systematic differences 
among the polar coordinates determined by different satellites. Details 
may be found in "Prediction of Polar Motion" by Y.S. Zhu, Dept, of Geodetic 
Science, 1981, in pr';.^ 

2.12 MERIT Short Campaign Data Analysis 

The BIN has already analysed the MERIT data using BIH Circular D as 
a common reference. Since new techniques are expected to obtain better 
accuracy than now achieved by the BIH, we think that a mutual comparison 
method may be more rational than using BIH Circular D as a common reference. 
In analyzing, we use both raw data and smoothed data. 

2.121 Using Raw Data 

Direct comparison . We compute the standard deviation of the polar 
coordinate differences of each two Analyzing Centers. The results are as 
follows: 



Q 

O 

X 

°Ay 

mean 

DMA - MEDOC 

0'.'040 

0'.'044 

0V042 

CNS(39) - CNS^o) 

C.031 

0.053 

0.042 

SAD - CNS^qj 

0.040 

0.036 

0.038 

CLA(0pt.) - SAD 

0.022 

0.025 

0.024 

CLA(0pt.) - DMA 

0.022 

0.025 

0.023 

CLA(0pt.) - UTX 

0.019 

0.026 

0.023 

DMA - SAD 

0.010 

0.011 

0.011 

UTX - SAD 

0.014 

0.012 

0.013 

UTX - DMA 

0.014 

0.010 

0.012 

SAD - CNS(3gj 

0.009 

0.014 

0.012 

UTX - CNS(39) 

0.010 

0.017 

0.014 
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From the above values we suggest that the standard errors of each 
Analyzing Center are approximately: O'.'O'i for MEDOC and for CNS(io)» 0*02 
for classical (optical), O’.'QOS - OVOl for SAO, DMA, UTX and CNS(3g). 

Relative polar motion comparison . Relative polar motion comparison 
is used to detect the possible systematic error (other than constant system- 
atic error— the difference of origin). Systematic error will be mostly 
reduced in relative polar motion, while it remains in direct comparison. 

The standard error of relative polar motion along with the standard 


error obtained in direct 

comparison 

are given 

below. 





Relative 



Direct 



a 

a 

a 

a . 

a. 
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AX 

Ay 

mean 

AX 

Ay 

mean 

DMA 

- MEDOC 0V039 

0'.'039 

0'.'039 

0'.'040 

0V044 

0'.'042 

SAO 

- UTX 0.010 

0.010 

0.010 

0.014 

0.012 

0.013 

DMA 

- UTX 0.014 

0.008 

0.011 

0.014 

0.010 

0.012 

DMA 

- SAO 0.010 

0.009 

0.010 

0.010 

0.011 

0.011 

Opt. 

- SAO 0.017 

0.022 

0.020 

0.022 

0.025 

0.024 

Opt. 

- DMA 0.024 

0.022 

0.023 

0.022 

0.025 

0.023 

Opt. 

- UTX 0.017 

0.024 

0.021 

0.019 

0.026 

0.024 

SAO 

- Circular D 0.005 

0.004 

0.005 

0.009 

0.009 

0.009 

UTX 

- Circular D 0.011 

0.010 

0.010 

0.011 

0.013 

0.012 


Since every is more or less larger than '’pgigt-jyg* ii^ means that 

there are systematic errors among the Analyzing Centers. The detailed 
properties of the systematic errors (say, the period, amplitude, etc.) 
could not be obtained because of the short time span of data available. 

2.122 Intercomparison After Smoothing 

In the initial stage direct comparison gave reasonably good evidence 
that systematic differences exist between the results of the earth rotation 
parameters obtained from different techniques. Later a smoothing process 
(by adjustment) was used by adopting a suitable mathematical model for inter 
comparison of pole positions obtained from different methods during the 
MERIT Short Campaign. The mathematical model used is given below. 
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X = ki + k2 cos A + ka sin A + k^ cos c + ks sin c 
y = ks - k2 sin A + ka cos A - K4 sin c + ks cos c 

where A = 2PI (MOD - 42413 )/365 
c = 2PI (MJD - 42413/435 
ki to ks are coefficients in seconds of arc 

The smoothed x, y coordinates of the pole were plotted in each case, and it 
was seen that curves depicting the motions were of different shapes and were 
at different locations with respect to the origin (CIO). These plots are 
Figs 1 to 7. 

As a consequence of the smoothing process, the values of coefficients 
kiandks obtained were plotted on a graph with their standard deviations 
(Fig. 8). ki and ks correspond to x and y coordinates of the origin of the 
curve, depicting the pole movements with respect to the CIO. This figure 
shows that systematic differences exist in the pole origin recovered as a 
result of these different techniques, which may be due to differences in 
mathematical models for computations and locations of observing stations, etc. 
Figs. 1 to 7 also give preliminary evidence that there is a systematic differ- 
ence existing in the values of annual and Chandler period amplitudes recovered 
by different techniques. 

The data available for the MERIT Short Campaign was not really sufficient 
to enable us to arrive at any conclusive evidence regarding the systematic 
differences and their estimation; and the noise level was also quite high. 
Therefore, we should wait for the intercomparison of data to be made avail- 
able through the MERIT Main Campaign. 
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2.2 VLBI Investigations - On the Time Delay Weight Matrix 
in VLBI Geoiletic Parameter Estimation 


2.21 Introduction 

A VLBI baseline consists of two antennas simultaneously observing the 
random signals emitted from an extragalactic radio source. The reception 
of the signal at one antenna is delayed in time fro.Ti its reception at the 
other antenna due to the difference in path length that the signal must 
travel. Denoting the time delay by t-j we have 

^ij ■ ■ ^i 

where t. (t-) is the time of arrival at station i (j). For a set of time 

‘ J 

delay observations about a triangle of stations we have tia, T 13 and T 23 . 
Assuming that we vere able to directly measure the times of arrival of the 
signal at each antenna (with equal precision) then it can easily be shown 
[Bock, 1980] that the unsealed variance-covariance matrix for that set of 
time delays is given by 

’l i -i 

i 1 i 

.-i i I 

Since the determinant of this matrix is zero, the three time delays in this 
ideal case are seen to be linearly dependent being related by 

■^12 + ^23 = 1^13 



Therefore, for each set of observations from a triangle of stations, one 
time delay would be redundant, and, furthermore, any two time delays would 
be correlated by a factor of i. However, in reality, the time delays are 
not estimated in this manner but are the result of a complex estimation 
procedure [Whitney, 1974; Whitney et al . , 1976] whose first step is the 
cross-correlation of the signals recorded on magnetic tape. In order to 
estimate, for instance, the pair xiz and T 13 , the signals from station 1 
are correlated with both the signals from station 2 and the signals from 
station 3; the station 1 signals, therefore, are involved in both cross- 
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correlations. The focus of this investigation is to assess the magnitude 
of the correlation between the cross-correlations on two baselines with a 
common station when all observations are performed simultaneously. If 
there is an indication of significant correlation (obviously the correlation 
will lie between 0 and i), the weight matrix for time delays used in the 
least squares estimation of geodetic parameters must be modified accordingly 
(in VLBI work today the weight matrix is assumed to be diagonal). See 
[Bock, 1980] for simulations on the effect of correlations on parameter 
precision estimates. The goal of centimeter accuracy in baseline determination 
for the detection of, for example, plate tectonic motion is the rationale 
behind this study, i.e., all possible modeling errors need to be investigated. 

This is a shortened and modified version of the final report to be 
presented at a later time. 

2.22 Group Delay Estimation and Statistics 

This very brief review of the measurement process is given through to 
the estimation of group delay (which is the estimate of the time delay so- 
called due to the fact that its estimation involves observations over a 
band of frequencies). The signal is collected at each antenna (as a voltage 
—see the discussion below on thi statistical model), passed through the 
receiver chain where it is corrupted by noise. The receiver chain includes a 
radio frequency amplifier which increases the power of the signal, followed 
by mixing the signal with a strong local oscillator signal, generated by 
a hydrogen maser, to convert the signal to an intermediate frequency. The 
resulting signal is again amplified and distributed at 28 narrow bands of 
2 mHz width [Rogers, 1979]. At this point the new video frequency signals 
in each channel are in analog form. The signals are clipped, sampled at 
the Nyquist rate, formatted and recorded on magnetic tape along with 
accurate timing codes. The clipping can be represented as a function 
which normalizes the amplitude of the signal voltage, X(t) such that 
[Thomas, 1976] 

F(X(t)) = *1 if |X(T)i > 0 

F(X(t)) = -1 if iX(T)i < 0 

The recorded signals are thus recorded in digital form. 
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The signal X 2 (t) received at the second station is equivalent to 
the signal X 2 (t) received at the first station but delayed in time by xi 2 
such that 

X 2 (t) = Xi(t - T12) 

The Fourier transform of X 2 (t) yields in the frequency domain 
xzioj) = / X2(t) e'^“^ dt = Xi(w) 

—CD 

where u is the radio frequency. 

It has been demonstrated by [Rogers, 1970] and [Whitney, 1974] th^^- 
the maximum likelihood estimate of (group) delay and delay rate (the time 
derivative of delay) can be obtained by maximizing the counter-rotated 
cross-spectrum (as defined in [Whitney, 1974]) summed over BT components 
where B is the channel bandwidth and T the total integration time, as 

max over Xi(u)) X*(w) ^^^t) + I 12 ] 

where T 12 , ti 2 » ii 2 are trial delay, delay rate and fringe phase based on the 
best a priori information available, and w is now the video frequency. The 
maximization above is equivalent to the cross-correlation of the signals re- 
corded at both stations with fringe rotation, that is, shifting one-bit stream, 
on the basis of a priori information, in order to align the two-bit streams to 
a point of near-maximum correlation, followed by multiplication of the signals. 
That the above expression represents a cross-correlation follows from the convo- 
lution theorem of Fourier analysis (see, for example, [Bath, 1974]) which states 
that if 

fi(t) t — *• Fi(ai) and F 2 (t) — »■ F 2 (u) 

then 

fl(t) a f2(t) « fi(a.) f2(ui) 

That is, convolution (symbolized by 0) in the time domain is equivalent to 
multiplication in the frequency domain. Convolution is defined by 

fl(t) a f2(t) = j fi(t) f2(t-T)dl = j"" fi(t-x) f2(T)dx 
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whereas cross-correlation is defined as 




fl(t-T)f2(t)dt 

. 


Thus convolution is the same as cross-correlation except that one of the 
time series should be taken in reverse order, i.e.. 


= fx(-T) 8 f2(T) 

where t in the above equations refers to the familiar time lag between two 
time series. Algorithms designed to approximate the cross-correlation, in 
practice, are given in [Whitney, 1974]. 

An expression for the normalized analog cross-correlation function 
assuming that both stations possess a square bandpass filter of width W is 
given in [Thomas, 1972a, b and 1973] as 



where is the model delay by which the signals are shifted in order to 
bring the data streams at both stations in close alignment. The expres- 
sion in brackets is the amplitude of the cross-correlation function for an 
extended source (a source that subtends a finite solid angle as opposed to 
a point source). The fringe visibility y accounts for the power lost due 
to self-interference of the extended source (y=1 for a point source); 

Wjj/W accounts for power lost due to imperfect bandpass alignment; the ratio 
of antenna temperatures T^ and T^ and receiver temperatures 
described in the next section account for signal-to-noise properties. The 
function indicates the accuracy with which the two signal streams have 
been aligned and peaks for at = t-T|^ = 0, t being the true time delay (includ 
ing atmospheric and instrumental delays as well as the purely geometric 
delay). The so-called fast fringes cos<>^ expresses the overall phase 
behavior of the cross-correlated signals. For digital signals, in the case 
when the correlation amplitude is small ('v- 0.001 - 0.1) as is the case in 
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VLBI work today, the digital cross-correlation function has a further 
loss of 2 /tt in correlation amplitude (see [Thomas, 1972b] also for the 
strong signal case). 

Once optimal values for delay and delay rate have been determined 
(for each frequency channel), the fringe phase (the phase difference) for 
the same frequency channel at stations 1 and 2 can be estimated by 
[Whitney, 1974] 

Im I Z Si 2 (w)j- 

tan $12 ? r 

Re I z 5 i2 ((*') ] 


where Si 2 (w) is the counter-rotated cross-sprectrum above using the optimal 
values of delay and delay rate for one channel. (Im is the imaginary part. 

Re the real part of the summations which are complex quantities.) 

In the bandwidth synthesis technique [Rogers, 1970, Whitney et al . , 1976], 
the fringe phase is determined over several frequency channels (spanning a total 
bandwidth of up to 400 mHz in the Mark III system) in order to improve the 
accuracy of group delay (and delay rate) estimation, and to resolve 2:: 
ambiguities in the fringe phase. The final estimate of group delay is the 
slope of a line passed in the least squares sense through the estimated 
phases of several frequency channels [Whitney, 1974]. The frequency channels 
are not generally uniformly spaced across the total bandwidth. Suppose we 
have observed a radio source from three stations, and that fringe phases 
have been estimated over n channels identical on each of the three baselines. 
Further suppose that the full variance-covariance matrix for the estimated 
phases over all channels is available. We shall assume that it takes on the 
following block-diagonal form 


T 



'(j)2 


0 



z 


n 


3n X 3n 


where 




0 ^ 0 

$12 

0 . 

'?12»‘f23 

E 

<>l 


4>i 3 

a 

4>i 3 3 


sym 


2 3 J 

Let us represent 






Pix P 12 

Pi 3 ' 


Pi = J:"' 


P 22 

sym 

P 2 3 
P 33 . 

; (P 
i 


block diagonal ) 


We have assumed that all fringe phases for one baseline (ij) have equal 
variances and are uncorrelated, though for each channel (observed simul- 
taneously by all three stations) the fringe phases are correlated. The 
group delays as mentioned above are estimated as the slope of phase versus 
frequency. We shall assume that the spacing of che frequency channels is 
large compared to the bandwidth of a single channel so that only the spacings 
between the frequency channels need to be considered [Whitney, 1974]. In 
order to determine the correlations between sets of time delays from three 
baselines we use the (linear) mathematical model of a line passed through 
the points refers to the video frequency of the channel, 

'*'ijk ff'inge phase on baseline ij, channel k. We form three sets of 
observation equations, one for each baseline 


♦ l2k 

+ 

Vl2k = 

Il2U>k 

+ 

bi2 

^3k 

+ 

Vl3k = 

TlJ^k 

+ 

bi3 

^23k 

+ 

V23k = 

I23^k 

+ 

b23 


The parameters include na, th, T 23 , the group delays to be estimated over 
all frequency channels and the b.., the line intercepts which for this dis- 
cussion can be regarded as nuisance parameters. The v^jj^ are observation 
residuals. We can write the above equations in matrix form as [Uotila, 1967] 

V » A X - L 


where 
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The least squares solution vector for the estimated parameters is given by 
X » (a'^PA)”^ A^PL 


In this discussion we are only interested in the variance-covariance matrix 
of the parameters (assume unsealed) 


“ (A^PA)‘i 


where the weight matrix P has been defined above. Performing the matrix 
multiplication we arrive at 
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ZuJ^Pll 

2cJ^Pi2 

“‘^|^Pl 2 

3 

2ii)|^Pl 3 


i^pll 

2 

^ P 12 

3 

2 Pl3 



Zw^p22 

Zu)|^P22 

2u^P2 3 

2 i*;|^p 2 3 

z, = 

X 



Ep22 

2oI|^P2 3 

J^P2 3 





Zu)^P2 3 

2oa|^P3 3 


sym 




SP 33 


The z. portion can be picked out after inversion including the time delay 
correlations. 

In the case of only two frequency channels uii, ^2 (say, for example, 
the two outermost channels (see discussion in [Thomas, 1973]), the group 

delay portion of the above z. matrix reduces to 
2 ^2 

k=l ^“k k*l 


j ^ a, , 

912 |^» 923 |^ 


Z. 

T 


1 

(a)2-wi)^ 





'k 


sym 



In the cases where the fringe phases from channel i for all three baselines 
are uncorrelated, i.e.. 


P 


1 1 





P^j = 0 (i j) 


the group delays are obviously also uncorrelated. In this case the variance of 

0 ^ , for example, can be determined from the first block of the now block 

■C 12 ’ 

diagonal z. matrix as 
X 
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*1 



in agreement with [Whitney, 1974] and [Thomas, 1976]. An expression for 
o’ is derived in [Whitney, 1974] (for T.<<To) 



where the notation has been explained previously. Furthermore, he defines 
the signal-to-noise ratio as 


SNR.. = 00 /51T 

where oo is the amplitude of the cross-correlation function described 
earlier (here given for an analog signal and a point source) and thus 
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2.23 The Nature of Signal and Noise In VLBI Observations 

As in all statistical detection problems, we are interested in the 
reception of a signal which bears information of value and noise which 
opposes cjr efforts to extract it. (This paragraph will draw from 
[Steinberg, 1963; and Kraus, 1966]). In radio astronomy both the signal 
and noise are of an equivalent statistical nature. Consider a resistance 
of R ohms at a temperature T (“K). The thermal agitation of the electrons 
in the resistor produces a voltage Og(t) at the terminals with zero mean, 

Ejn (t) * 0 and a flat spectrum Sn (ui) « 2kTR (k is Boltzman's constant). 
This so-called thermal or background noise can be shown experimentally to 
be a normal process, A resistance R is therefore equivalent to a source 
of voltage with zero mean and power, 2kTR in series with a noiseless 
resistor R. If such a ,oltage source is connected to an external load 
resistance R, , it will ieliver a certain power to it which will be maximum 
if the load resistence is equal to that of the source. The power delivered 
will be equal to kT. The voltages involved in VLBI antennas are of the 
same statistical nature as background noise. The same noise power is avail- 
able at the terminals of an antenna (kT). The antenna has a radiation 
resistance R and is exposed to a sky temperature, T. Therefore, we see that 
observations in radio astronomy consist of the measurement of an ante^ma 
temperature and the conversion of that temperature to a measure of power 
received from the source. The signals are corrupted by noise as they pass 
through the receiver components primarily at the first stage of amplifica- 
tion [Thomas, 1976]. This noise is basically background noise as described 
above. Thus, signal and ncise in VLBI observations have similar statistical 
properties, being modeled by the thermal noise of a resistance raised to a 
certain temperature. The receiver noise is expressed in terms of the 
radiation resistance of the antenna. It is the temperature to which the 
radiation resistance of the antenna would have to be raised in order to 
produce the same noise power as the complete receiver. The sum of antenna 
temperature, T^, and the receiver temperature, Tj^, is called the system 
temperature. Thus, the overall noise power is given by k(T^+T^). In order 
to appreciate the temperatures involved, for the Mark III wideband receiver 
and a 10 Jansky source (1 Jansky ■ 10 watts/m Hz) for the Westford 18 m 
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antenna, ■ 0?45 K and for the 4 m portable antenna, * 0?C3 K. The 
system temperature, on the other hand, is about 160® K [Rogers, 1979]. Even 
for large antennas observing strong extragalactic sources, the temperature 
wll 1 not usually exceed 1° K. Thus, antenna temperature Is about two orders 
of magnitude smaller than receiver temperature. The signal Is thus propor- 
tional to antenna temperature and the noise to receiver temperature. In 
order for the signal to be detected, measurements rmjst be Integrated over 
several minutes of time and the spanned bandwidth Increased (up to 400 mHz 
with the Mark III receiver). 

Considering the statistical nature of signal and noise in VLB! obser 
vations. It is reasonable to assign the following statistical model fjr 
output voltage of the receiver as [Rogers, 1970] 

X(t) « (T^)* s(t) + (T^)* n(t) 

where s(t) is the noise signal collected by the antenna from the radio 
source emission and n(t) is the receiver noise. The temperatures T^ and Ijj, 
therefore, serve as weighting factors. 

The signal and noise are complex quantities whose real and imaginary 
parts are assumed to be normally distributed (N (0,i)).1or a particular t, 
and independent [Whitney, 1974]. Therefore (and similarly for n(t)) 

E{s(t)} ■ E{s (t) + i s.(t)} ■ E{s_(t)) + iis.(t)} • 0 
r I r I 

var[s(t)] » E{|s(t)|^} « E{s‘(t) + st(t)} - E{s‘(t)} + E{s?(t)}» i+i*l 
cov[s(t) ,n(t)]» E{s(t) n*(t)} - E{s(t)} E{n*(t)} ■ Eisit) n*(t)} 

* Ets^(t) n^(t)} + E{s^.(t) n^(t)} + i[E{s^(t) n^(t)} 

- EiSj.(t) n.(t)}] • 0 

and the signal and noi-iC* are seen to be independent standard normal, N(0,1) 
random variables for 1<Aed t. Both s(t) and n(t) are assumed to be station- 
ary random processes (see [Swanson and Mathur, 1968]) so that their respec- 
tive autocorrelations are a function of time difference and independent of 
time origin. 

Since an observation Is integrated over a certain time interval T 
and a spanned bandwidth B, s(t) and n(t) can be represented over T by 
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a Fourier series expansion with BT Independent components. Furthermore, 
since the signal X(t) Is generally received from an extended source, we 
can assume that only the correlated part of the signal contributes to 
the antenna temperature and the retriainder contributes to the receiver 
temperature [Rogers, 1970]. 

The statistical model can just as well be written In the frequency 
domain as [Whitney, 1974] 

XU) • (T^)* sU) + (Tp)* n(uj) 

Both signal and noise are bandllmlted which means that their spectra equal 
zero outside certain regions of the frequency o [Papoulis, 1965]. The con- 
sequence of the above assumptions Is that the output voltage of the 
receiver X(t) at each station Is Itself a stationary, normal (with zero 
mean), bandllmlted stochastic process. 

2.24 Correlation Analysis 

We are now In a position to assess the degree of correlation between 
the cross-correlation of signals on any two baselines when the signals from 
one station are common to both processes. We shall Investigate a triangle 
of stations observing a radio source simultaneously(since this Is the basic 
unit for such an analysis) whose signals are represented statistically 
in the frequency domain as 

Xi(^) » (T^J* Si(w) + (TpJ* ni(^) 

X:(c.) - (T^J^ s:(^) + (T^J^ ni{^) 

X,(^) - (T^^)* S 3 (.) + (Tp^)* .i,(.) 

where 

s.(-) » si(.) 

s,(w) . sx(^) 

The notation is the same as in the previous section. In additioi., we assume 
that the noise introduced at each receiver is uncorreUtcd with that of 
another receiver, i.e., 

E (n.nj) « 0 (and E is.nj) • 0) (i^j) 
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we shell hetennine the degree of correlation between the cross-correlation 
on baseline 1-2 and 1-3. The rotated video cross-spectrum for baseline 
1-2 and baseline 1-3 can be written as 

SiiU) • Xif.;) X*^) e‘'“‘-‘"* 

Su(w) = Xi(w) 

which can be written in terms of signal and noise as 

+ Sx(ui) ni*(ui) + /TrJr7 

Ss 3 (w) - st(.) e**-t -'■(“) t 

+ /fTTr” si(^) 


where 


\ n ^ ^ „iu(Xl2+ 

02 (dj) - Haluj) e 

03 (til) — n 3 (tij) e 


Tizt) 

Tx3t) 


we shall assume as in [Whitney. 1974 ] that n](») has the same statistics as 

nzlixi) and likewise for n 3 '(^) and 03 ( 0 .). 

The above Si.(o) and Si,(.) are evaluated as the sum over a large 

number of independent short time intervals spanning an observation and 
refer to one component in the observed frequency spectrum. The maximiza- 
tion, as mentioned before, is performed by summation over BT independent 
points in the frequency spectrum. We can represent this summation as in 

[Whitney, 1974 ] by 


^X2 = zXi(bd) X*(i.) 

= ^ 17 ^ S*(ai) + -^rJaz 

+ ssCa'l n] (oi) + »^rJr7 


Eni(u)) s (w) 

* 

Eni(o)) 02 


and 
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= ^sU) s*(a.) + zni(.;) s*(..) 

+ i:s(u)) nj (id) + Zrii(a)) nj (id) 

where we have dropped the subscript from the signal terms. We are only 
concerned with the observations that are common to all three stations. It 
can be shown easily that all the terms in the expressions S12 and Si3 are 
independent of each other and the real and imaginary parts of each term 
are independent [Whitney, 1974]. Thus, ^12 and ^13 are vectors in the 
complex plane, and are the sums of independent random variables. 

The correlation between these two cross-correlation processes given 
by ^12 and ^13 can be defined as 

= cov[^i 2 . ^13] 

{var[^i2] var[^i3]}* 

which is also a vector in the complex plane since in general the covariance 
between two complex random variables (in our case, the sume of uncorrelated 
random variables) is complex but the variance is real [Jenkins and Watts, 1968]. 
For the covariance between ^12 and ^13, by definition [Papoulis, 1965], 

COv[^12,Ri 3] = E{[^12 - E{^X2)][^t3 ~ E{^!3}]> 

= E{^i2 fi3) - E{^i 2} ElSta} 

where the derivation is given in the final report. To derive an expression 
for the variances we can use the property that the four terms of both ^12 
and ^13 are independent and therefore the variances of each term can be 
summed. By definition [Papoulis, 1965] 

^12 
and 

!tx3 


E{|^i 2 - El^izll’} 


E{1^13 - E{^i3 }i"} 

^T,^)(T,^ .T,^) 
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where the derivations are given in the final report. Then, for the correlation 
vector 




■EE 


(Ta," Tr.1 




The magnitude of 


Ipa ^ 1 


5E 




* 


is a measure of the degree of correlation between the two cross-correlation 

processes, represented by ^12 and ^13. (The same result has been derived 

independently by Herring [private communication].) The argun«nt of 0 * * 

•'12 ,f'i 3 

_ Im p 

arg Pjf jf = tan"^ — - = <j>i2 

Ri2,Ri 3 Re ^ 


Pi 3 


expresses the effective phase difference between the two processes-- 
4112 - 4 >i 3 is the angle between ^32 and ^13 in the complex plane. The 
magnitude of the correlation vector is defined in general by [Jenkins and 
Watts, 1968 ] as the coherence between two complex random variables, Xi and 
X2. by 

|cov [Xi, X*]l 

ki 2 = 

{var [Xi] var [X2]}* 


The correlation vector is also referred to in the literature as the complex 
correlation factor used, for example, in optics for the study of light from 
an extended, incoherent quasi-monochromatic source [Born and Wolf, 1964 ]. 

Extending the above results, the following Hermitian matrix (similar 
to the coherency matrix in optics) for the cross-correlation of signals in 
a VLBI triangle can be defined as 


J 


var [^32] 

COV [^13a^12] 

cov [^23.^12] 


COV [^12»^13] 
var [^33] 
cov [i? 23 .i^l 3 ] 


cov [^32,^23] 
cov [i^i 3.^23] 
var [I23] 
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'Sjsa 




g*i 

g-l ('J’12“<}>23) 




^Sx's, 

->3 "1 "2 


^S2*^aJa3 


^Sa^Ai'^Aa 


'S2^S3 


gi ('J>12”'!>2 3 ) 
gi(<f'l3“<^2 3) 


where T- = T. + T„ , the system temperature for station j. 

J 

The corresponding coherency (correlation) coefficient matrix is 
given by 


1 


^^^ 12 ,^ 23 ^ 




1 


= 

• IW. 



1 


1 . 


where we can pick out the coherence between any two cross-correlations 
from the off-diagonal elements. 

It is interesting to note that the coherence between the cross- 
correlations for any two baselines (with a common station) is just the 

amplitude (denoted here by IpL. ) of the normalized cross-correlation 

^ j k 

function defined earlier (or as is defined in [Whitney, 1974], the true 
correlation coefficient of the third basel ine) ,i .e. , 





for the case of observations to a point source with the signals recorded in 
analog form. 

In order to assess the magnitude of the off-diagonal elements of 
matrix C, i.e., the degree of coherence (correlation) between cross- 
correlations on two baselines with a common station, we recall from the 
earlier discussion that in VLBI work today This indicates that 

the C matrix is very nearly the identity matrix, i.e., the cross-correlation 
pairs are virtually uncorrelated (|ph 0.001 for the example of the previous 
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section). It follows that the fringe phases are uncorrelated and consequently 
the group delays (as well as the delay rates). In fact, the correlations can 
be expected to be even smaller considering that we are observing extended 
sources so that the correlation would have to be multiplied by a fringe 
visibility factor (remember 0<y<1). In addition, we know that digitizing 
of the signals reduces the amplitude of the cross correlation function by 
a factor of Zh and (in the case of weak signals) can be expected to affect 
the cor'* Jations. Finally, the theoretical cross-correlation procedure 
as described earlier is only approximated in practice, which will also tend 
to reduce the degree of correlation (see [Rogers et al., 1979]). 

On the other extreme, in the limiting case where 

C = 

as could be envisioned for satellite interferometry where a radio antenna 

(e.g., SERIES) would observe an artificial satellite transmitting a much 

stronger signal than extragalactic radio signals. Furthermore, the satellite 

could be considered a point source and digitization would not be expected to 

decrease the cross-correlation function amplitude (see [Thomas, 1972b]). It 

then could be inferred (this will be shown in the final report) that the 

fringe phases are correlated in the same manner. Since the C matrix in 

this case is singular, this implies that only two cross-correlations in a 

triangle are independent and likewise for the fringe phase. The covariances 

of the matrices defined earlier could then be formed by 
^i 

where p, the correlation between fringe phases will approach i. It is easy 
to see (at least in the two-channel case) that the correlations between 
group delays would approach that of the hypothetical situation described 
in the introduction and could be computed in general by the portion of 
the matrix developed in Section 2.2.2. 


1 i r 

1 i 
sym 1. 
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2.3 Utilization of Range Difference Observations in Geodynamics 


2.31 Utilization of Lageos Laser Range Differences 
Introduction 

The objectives of this project as originally stated in the proposal 
are the development and implementation of the technique of range differencing 
with Lageos in order to obtain more accurate estimates of baseline lengths 
and polar motion variations. It is expected, and the simulations done to date 
confirm this, that by means of differencing quasi -simultaneous range observa- 
tions a great deal of orbital biases can be eliminated resulting in an esti- 
mate which is virtually bias free. In the past few months attempts were 
made to use real data in order to assert the conclusions made on the basis 
of the simulations. Partly due to the oversimplified models used in our 
software (GE0SPP80) and partly due to the very poor geometry and distribution 
of the data, these attempts fell short of providing any conclusive results. 

It was realized, however, that even with the best possible geometry and dis- 
tribution of the observations, certain physical phenomena would have to be 
included in the model. This being the case, a complete revision of our 
software was undertaken, which resulted in a new version (GE0SPP81) cur- 
rently undergoing tests. Given in the following is a brief summary of the 
various models implemented in the new version. 

2.311 Data Preprocessing 

As mentioned in previous reports, it is quite improbable, if not 
impossible, to obtain exactly simultaneous laser observations from two 
ground stations to a satellite. Modern lasers, however, have high repeti- 
tion rates and given fair weather conditions and accurate predictions, a 
sequence of ranges at a rate of about 1 pps can be easily achieved. Since 
an average Lageos pass lasts about half an hour, this implies a large amount 
of data. The high altitude of the target makes it possible to observe it 
from several stations simultaneously, even if the stations' separation is 
of continental extent. The current procedure to obtain the quasi-simultaneous 
ranges from data sets such as described above is to determine the "overlap" 
observations for the station pairs involved, isolate these observations and 
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determine which of the two stations in each pair has a denser data distribu- 
tion. Once this is determined, the data of the station with the most obser- 
vations are fed into a cubic-spline interpolator and ranges are obtained for 
each of the data points in the alternate station's record. These ranges 
then are differenced to produce the range-difference data for GE0SPP81. The 
data used in this procedure have already been corrected for systematic errors 
(Tape-90 data in GEODYN terminology). Figs. 9 and 10 depict the data record 
for two passes of Lageos. The bars indicate the epochs when the actual ob- 
servation occurred and the curve which joins their centers is the spline fit 
to these data. The stations marked with {*) are the ones for which we used 
the original observations in the formation of the range differences. 

2.312 Reference Frames 

A number of reference frames are involved depending on the phenomena 
being treated in the program. The integration of the orbit is performed in 
an "inertial" frame defined by the mean equator and equinox at some specified 
date (reference epoch) and the atomic time scale (TAI). Perturbations are 
computed in "convenience" frames, which will be explained later, and then 
rotated into this frame in order to evaluate the equations of motion and 
the variational equations of state. 

The position of the observing stations is defined in a geocentric 
earth-fixed system which is materialized by a mean pole (e.g., CIO) and a 
mean Greenwich meridian (e.g., BIN). The observations are tagged with UTC 
epochs. The link between this system and the inertial system is provided 
by the Greenwich hour angle of the true vernal equinox and the coordinates 
of the true celestial pole with respect to the mean pole used. The link 
between the inertial frame for the integration of the orbit and the true 
of date frame at the observational epochs is provided by the precession 
and nutation theories of Newcomb [ESAENA, 1974] and Wahr [1979] respectively. 
The time scales involved, UTl and UTC are related to each other using the 
difference a(UT 1 - UTC) as obtained from JPL [Fliegel, 1981]. 

2.313 Gravitational Perturbations 

The perturbations due to the Earth's gravitational field were discussed 
in the previous report. Although that presentation is basically valid. 
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thire is a subtle point in the coordinate system transformations involved 
which deserves further explanation, especially in view of the existing con- 
fusion In the geodetic community. This will be dealt with in the next 
section. The new version of the program includes point mass accelerations 
due to the attraction of the moon, the sun, and the planets Venus, Mars, 
Jupiter and Saturn specified at option by the user. The planetary ephem- 
eris, as well as the lunar and solar position vectors, are obtained by 
interpolation from the JPL tape containing the DE96 Solar, Lunar and Plane- 
tary Ephemeris. 


2.314 Computation of the Non-Spherical Effects from the Terrestrial 
Gravitation 


The gravitational potential is usually expressed as a series expansion 
in spherical harmonics. With the 0‘'^-order term (point mass effect) removed, 
we can write the non-spherical part of it as: 

«e ^ 


V(r,^,X) » Y 





nm’ 


( 1 ) 


where r.iii.x are the spherical coordinates of the point of evaluation. From 
the purely mathematical point of view the choice or definition of the (r,4i,.\) 
system is irrelevant. The function V, the potential, which is approximated 
by the series is invariant with respect to any similarity transformation of 
the underlying coordinate system; the value of V at P(r,<|),A), Vp will be the 
same whether P is defined in the system or, say, the (r‘, 4 ) 

system. The series, however, in (r',(j(‘,A') will have different constants, 
harmonics, in this case. It is thus obvious that the similarity transforma- 
tion between (r,^,\) and. (r‘,i?',x') propagates as similarity transformation 

between (C„, S„) and (C' , S'). 

nm nm nm nm 

So far, the above is just a restatement of already known facts. For 
instance, Kleusberg [1980] has developed the (rather cumbersome) formulae 
to obtain the primed harmonics from the original, given the transformation 
parameters. It is the last statement in the previous paragraph that causes 
the problems: The fact that a change in the coefficients cannot be attributed 
to an actual change of the coefficient or a coordinate system cnange, unless 
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it is known a priori which of the two happened. In a sense this is similar 
to "the principle of equivalence" in general relativity, whereby gravitational 
and inertial forces are inseparable. Considering now that in a dynamical 
solution the coordinate system definition is provided by the satellite 
dynamics (short of a longitude definition), one realizes how important it 
is clearly define a priori the system in which the harmonics of the 
series in (1) are referenced. From the geodetic point of view the determina- 
tion of these hamronics, which is accompanied by a determination of station 
positi -.s for a large number of globally distributed stations, provides by 
itself definition of the underlying coordinate system via the coordinates 
of these stations. When using such a set of harmonics in a different problem, 
(e.g., in a partial solution for new station locations), we have two options 
as to maintaining the internal consistency of the solution: (a) either 

include observations from the fundamental stations which are constrained to 
the original positions (very impractical if not impossible), or (b) know 
a priori the relation of the reference frame for the geopotential to any 
other frame involved in the solution. The second option, which is practically 
the only one available to us, will be examined in the following, in connection 
with the evaluation of the equations of motion and the variational equations 
of state. 

Since the geopotential is given in an earth-fixed system (r, 4 i,x), it 

D 

is simpler to compute the required derivatives in that system and rotate 
the results into the inertial frame of integration: 


” I 

%S 

- I 

where 

''b 

M 

S 

N 

P 


= (SNP)^ M Vb = CP^ s"^ M] (2) 

inertial acceleration due to non-spherical effects 

= vV(r,iii,x), body-fixed gradient of the potential in (r,<ii,x)|^ 

the Jacobian which rotates from (r,<j>,x) to the (x,y,z) system 
earth orientation rotation for GAST, ^p>yp 

nutation matrix for the celestial pole 
precession matrix for the celestial pole 


The N and P rotations are straightforward and need not concern us anymore. 
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The matrix S rotates the true of date system into the mean earth-fixed 
system (x,y,z)|^ in which Xp and yp denote the polar coordinates of the 
former. The transpose simply reverses the sense of the transformation. 
This then implies that the vector M is in the mean earth-fixed system 
(x,y,z)jj, which in turn implies that the matrix M rotates into {x.y.z)^^, 
a vector defined in the ( r, «}>, A )b system, consistent with the expansion of 
V in ( 1 ). If this system happens to be the mean earth-fixed system, then 
M takes the usual form: 


M 


-xz 


r^*^+y" 


-Jl. 




--.y 

x^+y^ 

X 

x^+y^ 

0 


X 

r 

I 

r 

z 

r 


(3) 


where (x,y,z) are the mean earth-fixed coordinates of the evaluation point. 

If on the other hand the spherical harmonic expansion in ( 1 ) is defined in 
a system (r' ,(}>' ,x' ) , then 

Vjj = M-V- (4) 


which implies the following: 

(a) The potential gradient (r',4>',x') will have to be evaluated using 
the satellite spherical coordinates in (r' ,41 ' ,x' ). rather than (r,4i,X).. 

(b) In order to use ( 2 ), will have to be rotated as in ( 4 ), to become 
consistent with the rest of the rotations. 

In the general case, therefore, where one accepts that the geopotential 
refers to a system (r',$',x‘), other than the one implied by the elements of 
the S matnx, equation ( 2 ) must be rewritten as follows: 


" I 


[p’’’ M M'] 


( 5 ) 


where M' is the Jacobian between (r,<ji,x) and (r',4',x’) 
is the gradient of V in (r',4',A') 
and the other matrices are as in {2} 

This problem of inconsistency becomes particularly important in the 
case when one estimates polar motion and earth rotation and then compares 
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these results to those obtained by other agencies. Even if the same 
data were used for all solutions, the use of different geopotentials 
withouc use of equation (5) will result in systematic differences in the 
results. 

The above comments apply in the same way in the computation of the 
non-spherical effects contributions in the variational equations: 


■ . r 
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In the computation of — I 


3r. 


, the follov/ing terms will 


NS 


computed: 




'_3i 

-n- 

» 




have to be 




( 8 ) 


where r, 4 >,A are the coordinates' arguments of the geopotential V. If these 
are not consistent with the (x,y,z) system in which r^ is referenced, but 
rather with (x',yiz‘) such that 


(x',y',z*) = L(x,y,z) (9) 

then the terms in (8) must be replaced by 


r 3r 1 


r3r ^ 
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and similarly 
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( 10 ) 



In a more compact form, since 
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using (9) and (10), we can write (11) in the most general form as 




3(r,({.,,\) 

1 r 

3r‘ 

u b J 

U(r,^,X)^^ i 

[/Vns 



3r 
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and denoting by 




equation (6) becomes in the general case 
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(13) 


or 
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(14) 


This discussion would hardly be complete if one did not identify where 
the real problem lies: Why geopotential models obtaineo from similar types 

and amounts of data differ in the coordinate system definition. The obvious 
reason would be that the stations which participate in the observations and 
define (with their estimated coordinates) the reference frame in each case 
are different. This is true, but in practice it is a large number of 
globally distributed stations which participate in such campaigns and one 
would think that on the average the geometry would have little effect, 
certainly not introducing biases with respect to certain regional sub- 
networks. The problem seems to lie in the neglect of certain physical phe- 
nomena, affecting the harmonics themselves, in the process of their estimation. 
It is a well known fact [Heiskanen and Moritz, 1967-, Nagel, 1976] that for 
the low-degree harmonics (up to n=2) we can easily associate with them 
certain properties of the estimated model. The pair (C,^, is of 
partivcular interest to us since their values are directly proportional to 
the alignment of the reference system with the principal axis of maximum 
moment of inertia. Limbeck first called attention to the Implications of 
this otherwise innocent looking property, in [Lambeck, 1971], where he 
stated that 
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u* = and V* S - ^ (15) 

t20 1-20 

u* and V* being the coordinates of the axis of figure with respect to the 
third axis of the reference coordinate system (the axes u*,v* are defined 
in the sane sense as the x ,x for the pole). A proof for (15) can be 
found in [Nagel, 1976]. This interesting property becomes a real problem 
in practice, due to the fact that the axis of figure is time variant. In 
[McClure, 1973; Nagel, 1976; Leick, 1978; and Moritz, 1979], one will find 
detailed descriptions of the motions of this axis. It suffices here to 
mention that there is a free motion with an amplitude of 2m and a period 
equal to that of the Chandlerian wobble ('^- 430 days) and a forced diurnal 
motion with an amplitude equal to 60 mi The elastic-earth-model motions 
for the various axes involved are depicted in Fig. 11, taken from [Nagel, 

1976]. Figs. 12 and 13 [Nagel, 1976] show the diurnal variation of C 21 and 
S 21 due to the diurnal motion of the axis of figure. Ir, Fig. 14 [ibid.], 
the diurnal variations have been filtered out and the effects of the free 
motion on C 21 and S 21 are shown, based on the BIH polar motion over 1968.0 
- 1973.0 period (from which the motion of the axis of figure has been deduced). 
It is obvious from these graphs that if one is trying to clearly define the 
coordinate system to which the estimated harmonics refer, one cannot afford 
to neglect such effects. Since the observations used in such major solutions 
span quite a long period of time (over five years in most cases), it is 
quite impossible to keep the coordinate system definition intact while at 
the same time solving only for average values of C 21 and $ 21 . As it was 
shown in this discussion, an improper coordinate system transformation in 
a later solution will introduce biases in the orbit and in the estimated 
parameters based on that orbit (e.g., polar motion). Based on the above, it 
seems that the only procedure which will produce confusion-free results is 
the one in which the timelike variations of C 21 and S 21 have been considered 
in the determination of the geopotential harmonics and the solution is 
performed for a prespecified epoch. In this sense, the orientation of the 
coordinate system in which the geopotential is referenced will be determined 
by the adopted model for the axis of figure motion, while the geocentricity 
of the system will still be maintained by the suppression of the first- 
degree harmonics. The motions of the axis of figure do not affect only 
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C;:i and S 21 , but the effects on the rest of the harmonics are well below 
the order of 10'^**, and for the time being they can be safely neglected. 


2.315 Perturbations from the Solid Earth Tides 

The model for the tidal potential is a rather simplified one account 
ing for effects to degree n=2 only: 

U, = p3 R* p2(cos4i) (16) 

I 2 ® 


Assuming that for an elastic earth k 2 s 0.29, the disturbing potential due 
to a body b at a satellite position r will be: 


Uoir) = 


kj GMj. 
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[3(ujj. u^)‘ 
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(17) 


where u. and u denote unit vectors in the directions of the disturbing body 
D r 

and the satellite respectively. The acceleration sensed by the satellite 
is = V Up(r), or from (17): 
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b(Ujj* u^)^] u^ + 2(uj^»u^) Ujj 


(18) 


The contribution in the variational equations is the following: 



10(ujj. u^)tu^ uj + Ujj u]^]| (19) 

and 



L3r J 

The above formulae must be summed over the bodies contributing in the 
tides. In 6E0SPP81, the lunar and solar effects are the only ones considered. 
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2.316 Solar Radiation Pressure Perturbations 


The model implemented in the program is a rather simple one, consider- 
ing only mean solar flux effects. The acceleration is given as 





V 


R.. 


VS 


m 




( 20 ) 


where 

V = shadow factor » 1- sunlight, or 0- shadow 

S = mean solar flux at 1 AU 

c = speed of light in vacuum 

Cp = radiation reflectivity coefficient for the satellite 
A = cross-sectional area perpendicular to direction of incidence 
m = satellite mass 

*^SUN ~ distance of earth-sun centers of mass (= 1 AU) 

7 = satellite position vector 

F 

s = sun's position vector 


From simple geometric considerations 


V = 1 if r»R > 0, or if r*l? < 0 but 
s s 
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From (20) the contribution to the variational equations is obtained as 
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Sections 5 and f are explained in detail in [Cappellari et al . , 1976]. 
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2.32 Doppler Experiments 
Introduction 


The Idea of using ranges derived through simultaneous Doppler obser- 
vations for a geometric solution has been fully explained in the previous 
semiannual report. It should be useful though to note that the data used 
for this experiment is the EDOC-2 campaign data set. There has been some 
progress in a small experiment involving observations from six stations 
which is the minimum for a geometric solution. However, an effort has 
just begun to perform the same solution with nwre stations involved and 
many more common observations included. Progress and problems arising 
during this effort will be included in this report. 

Additionally, there has been another Doppler experiment which started 
a few months ago. This is also a geometric range solution, this time using 
a different set of simultaneously obtained Doppler data. This data set is 
the Victoriaville data sent to us from the Ministere de I'Energie et des 
Resources of Canada. At the end of this chapter, p -ogress on this effort 
will be reported. 

2.321 Geometric Solution Using Ranges Derived Through Simultaneous Doppler 
Observations 

The Ranges 

The Doppler derived ranges have been evaluated according to the theory 
described in the previous semiannual report. These ranges which are the 
observations in our adjustment should be very close to the corresponding 
geometrically derived ranges. The geometrically derived ranges are those 
evaluated through the formula 

r = {(Xg - + (Yg - + (Zg - Zg)2)* 

where Xg,Yg,Zg are the earth-fixed Cartesian coordinates referring to the 
ground stations. Xg,Yg,Zg are the satellite's earth-fixed Cartesian coordi- 
nates derived through the broadcast ephemeris. 

Indeed, the Doppler derived and the geometrically derived ranges have 
been compared through graphs and their differences were not really significant 
especially in that part of the pass where the range decreases (the elevation 
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angle increases). AFter the satellite reaches the highest elevation angle, 
the aforementioned difference becomes larger (sometimes 5 km) but still with 
such differences the geometrically derived ranges can be used as approximate 
ranges in the adjustment. 

The Geometric Range Adjustment 

The mathematical model used for the geometric range adjustment is the 
following: 

r • [(Xj - Xg)2 + (Yj - Yg)2 + (Z5 - Zg)^]^ + a^ + f(p) (1) 

where a^ is a constant bias including the error to the initial range plus 
the error due to clock offset among the stations; f(p) is a function of any 
kind of parameters affecting the observations (Doppler ranges). We will 
call the + f(p) part of the mathematical model the residual model. We 
will call preliminary adjustment the part of the adjustment where we keep 
the station coordinates fixed and we solve for the satellite coordinates and 
the coefficients of the residual model. 

The next step was to come up with a residual model which yields conver- 
gence and the smallest after-adjustment residuals possible. The method used 
was the following: First we performed a preliminary adjustment without using 

any kind of residual model. Our mathematical model was 

r » [(X3 - Xg)2 + (Y3 - Yg)2 + (Z3 - Zg)2]* (2) 

Then we tried to study the behavior of the residuals from this adjust- 
ment. These residuals have been plotted and the plots have been compared 
to each other for all the stations and the passes. 

The fact is that the part of the passes used for the adjustment is 
really small (about 10 minutes) to give a good estimation of the behavior 
of the residuals, but in all the cases we dealt with they had more or less 
the following pattern: Relatively small in the beginning, increasing posi- 

tively or negatively with time, showing a second degree or sometimes a 
third-degree polynomial behavior. The uncertainty in choosing one residual 
model led us to try using different models and study the rapidity of conver- 
gence and the magnitude as well as the behavior of the new residuals. The 
residual models tested were the following: 
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(a) V • a^ ♦ b^t (d) , . a^ ♦ b^t ♦ c^r 

(b) V . a^ ♦ b^t t (e) V a^, * b^t » ♦ d^r 

(c) V ■ + b^t + Cq cosec e 


where v are the range residuals derived using the mathematical model described 
in eq. ( 2 ), t is the time of the particular observation, e is the corresponding 
elevation angle, and r is the range itself corresponding to the particular 
observation. 

From all these models, (b) has been proven to yeild the smallest post- 
preliminary adjustment residuals as well as to be the fastest converging model. 
The linear model (a) yields post-preliminary adjustment residuals fairly large 
(to the order of 100 m) and having sinusoidal behavior. The residual model (c) 
yields an oscillating solution in terms of determining the position of the satel- 
lite as well as in determininc a , b , and c residual model coefficients. 

* 0 0 0 

Finally, using models (d) and (e) we did not have convergence in our 
preliminary adjustment. Therefore, the mathematical model we decided to 
adopt is the following: 


r . [(X 5 - X/ * (Yj - x/ * (Zj - Zq) 2 ]‘ . t b„t + 


The oost-prel iminary adjustment residuals we get using this model are of the 
order of 0.5 - 1 m during the biggest part of the pass, but they systematically 
become large (about 10 m) towards the end of the pass. Presently we are trying 
to explain somehow the systematic behaviorwhich has not been eliminated 
using all the previously mentioned residual models. 

It has been mentioned in the previous semiannual report that 19 passes 
have been selected to be included in the final adjustment. Fig. 15 shows 
how the projections of these common passes are located with respect to six 
observing stations. Some of the passes are not close to the- observing stations. 
Their elevation angles are small and in the preliminary adjustment they con- 
verge very slowly in terms of satellite coordinates yielding very large resid- 
uals. There is one pass for which the common part starts right above station 
LGA in Spain. For this pass due to loss of information we have no conver- 
gence in the preliminary adjustment. 

After all these investigations we came up with nine good passes con- 
sisting of 5694 common observations totally. This is probably not a good 
number for a strong geometric solution, but the experiment will be performed 
with those few observations and depending on the results we get we will add 
more passes in the future. 
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2.322 The Extended Geometric Range Adjustment 

A new attempt began a few months ago in the Department of Geodetic 
Science to perform a geometric range adjustment covering a large area in 
Europe. The data set used for this purpose is also the EDOC-2 campaign 
data set. Performing such a large geometric range adjustment we will be 
able to compare the results derived from various dynamic solutions performed 
by several European agencies. For this purpose all the stations equipped with 
JMR and CMA receivers from which data is available, 23 stations totally, have 
been selected initially to participate in this experiment. 

The Data Editing 

The first step in data editing was to write a program which takes as 
input all the tapes as sent from Europe and perfonns the following operations: 

1. Identifies the passes 

2. Rejects bad observations. This has as a result the break of the pass 

into two or more parts. We characterize as bad observations: 

a) those which have a recorded time interval which is not 4.6 or 4.9 s 

b) every 25th observation which does not have a 4.9 s time interval 

c) every 26th observation which does not start with an even 2 minute 
mark 

d) observations for which the recorded Doppler counts are not within 
the possible limits 

The output of this program is the information for each observation needed for 
the computation of ranges. The output also includes information for each 
individual pass containing the meteorological data for the tropospheric 
refraction. 

The next step was to find the common observations. For this purpose, 
a new program has been written which gives the parts of the passes which 
have been simultaneously observed from at least six stations. Using this 
program we came up with more than 400 common passes with length varying 
between 2 and 10 minutes. All the passes have also been plotted with 
respect to time and station and from this diagram it is very easy to identify 
their common parts. 
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This experiment will continue in the same pattern as the small one 
including the six stations in Spain, and we hope it will give a complete 
image of the advantages and the disadvantages of the new method of simul- 
taneous Doppler ranging. 

2.323 The Victoriaville Data Set 

The last experiment involving a geometric I'ange solution is one 
performed with the Victoriaville data set including Doppler observations 
from eight stations. These observations took place between October 2-4, 
1979. These data have been used by Sheltech, Ltd. under contract with 
Quebec Lands and Forests for the purpose of evaluating the accuracy of 
Doppler survey methods for establishing network control. The participating 
stations' locations are given in the following table; 




X 

h 

1 

45”45'26'.'0040 

72^16 '38^797 

218.29 m 

2 

45 37.6413 

11 39.2460 

172.60 

3 

52 16.2391 

15 43.2894 

104.02 

4 

57 55.3552 

02 58.1076 

158.33 

5 

55 34.8419 

07 52.7973 

112.82 

6 

56 53.5586 

16 09.4327 

101.60 

7 

51 01.6013 

04 59.7844 

132.27 

8 

49 52.5628 

00 58.5315 

124.22 


The observed satellites were 30190, 30140, 30130, 30120 and 30200. 

The Victoriaville data have been sent to us on one tape including the 
Doppler counts, timing information, broadcast ephemeris information and 
signal strength information. All these data were already majority voted. 

From the broadcast ephemeris information of the tape, using H. White's 
program, the observed satellites' state vectors have been derived. For two 
of the observed satellites, specifically 30140 and 30190, the precise ephem- 
eris was available and therefore a comparison of the state vectors was possi- 
ble. Tne state vectors for these two satellites were very close which fact 
assured us of the efficiency of H. White's program. Meteorological data 
for this data set has been sent separately. These meteorological data have 
been recorded at the locations of the stations every ten minutes and include 
temperature, relative humidity and barometric pressure. 
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The data as they have been sent to us have been used a: Input in the 
program JMR OOP-A provided by N6S and have been transformed to a format 
useful for further processing. 

2.324 Doppler Intercomparison Experiment 

Work is still continuing on the programs which will be used to reduce 
the data from the OSU Doppler Intercomparison Experiment (which was described 
in the Fourth Semiannual Report). Serious problems were encountered whici) 
prevented the use of the CALSAGA program for the final data reduction. There- 
fore the program system GEODOP has been obtained fromthe Canadian Geodetic 
Survey for this purpose. Work is underway to convert this CDC Fortran 
program to IBM Fortran for use on our system and to allow input of our 
intercomparison data. 

Problems with GALSAGA 

As explained in the last semiannual report, the GA’.SAGA (or SAGA) and 
SAMVAP [Brown, 1973] programs and the preprocessing programs written here 
(SUPASS and SA6SET) were to have been used to process the intercomparison 
data. However, after numerous test runs with some JMR-IA data, it was 
apparent that the GALSAGA program would not converge to a reasonable solu- 
tion, These test runs included the use of; a) up to 40 balanced satellite 
passes, b) various program constraints; c) direct and SAMVAP created precise 
ephemeris state vectors to establish the orbit, d) various weights for the 
station, state vector, and program constraints. 

The best solutions obtained during these tests resulted in final 
station positions which were incorrect by several hundred meters, unreason- 
able values for orbital and other parameters, and values for v"^PV of the 
adjustment which were far too high, even after up to ten iterations of the 
solution. There seem to be two possible explanations for these poor solu- 
tions. The first is based on the observation that GALSAuA was intended to 
be used for range observations, even though it does contain input provisions 
for Geoceiver data. It is slightly possible that some error had been made 
in the programming here to allow it to use JMR (and other) data, although 
this is not likely since it has been used successfully here for multi- 
stations solutions of JMR data [Arur, 1977], with only slight modification 
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since then. The second explanation is far more likely, being that GALSAQA 
will only operate well for multi-station solutions, since it is after all a 
short-arc program. A single station solution might not converge unless all 
coordinates (station and satellite) and paranseters are very highly constrained. 
The fact that there appears to be no specific references to single station 
GALSAGA solutions supports this conclusion. 

However, some further use of the GALSAGA program is still being 
considered. It is possible that another attempt at a single station 
solution may be made (with different data). If a single station solution 
can be obtained GALSAGA solutions could then be compared with GEODOP solu- 
tion. More importantly, GALSAGA may be used to process multi -station solu- 
tions for the Columbus-Ottawa baseline for which we have data. (There 
presently is no reason to suspect that the OSU version of GALSAGA will not 
work with two or more station solutions.) Therefore some work is being done 
to modify the SUPASS and SAGSET programs to; 

a) Read majority voted data in a PREDOP compatible format, to allow solutions 
with either GALSAGA or GEODOP with the same data sets, 

b) Use "30 second" Doppler counts instead of "4.6 second" counts, again to 
maintain compatibility with GEODOP and to reduce the computational expense 
in running GALSAGA, 

c) allow SAGSET to input multi-station solutions to GALSAGA. It should be 
noted, however, that this work will be of low priority, especially if our 
GEODOP version will allow us to adequately answer the original problems (as 
posed in the Fourth Semiannual Report. 

GEODOP Program Sytem 

Since solutions could not be obtained with the GALSAGA program, Mr. uan 
Kouba of the Canadian Geodetic Survey (CGS) was contacted to determine the 
availability of his GEODOP program and assoicated subprograms (such as PREDOP, 
PREPAR, MERGE, NWLFIT, etc.), herein called the GEODOP program system [Kouba, 
1976]. Through Mr. Kouba, the CGS provided us with a tape copy of the latest 
versions of the programs in CDC Fortran, along with the program manuals and 
updates to them. Since it would be necessary to fully understand the 
programs' operation, and to convert them to IBM Fortran, Brent Archinal of 
OSU traveled to Ottawa to confer with Mr. Kouba and others about the program 
system and to run some of the OSU data. These computer runs would then be 
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available for comparison with runs here to check our IBM version of the 
program. 

The conversion to the IBM Fortran is now well underway, with the 
PREDOP and GEOOOP programs now in operating condition here, except for a 
few problems still being worked on. The program system presently accepts 
CMA-751 majority voted data, with subroutines available (but still being 
tested) to use MX-1502 raw oata, and JMR-IA majority voted data. In the 
near future, we expect to: a) Finish the PREDOP and GEOOOP conversion, and 
testing of the input subroutines, b) test our conversions of the other 
utility programs in the system (MERGE, NWLFIT, etc.), c) Design the JCL 
language to operate the programs with the large amount of data (files) which 
will be used. 

After these conversions are tested and complete, the Columbus and 
Columbus-Ottawa Intercnmoarison data will be processed as explained in 
the Fourth Semi annual Report, the eventual goal being the determination of 
the most accurate instrument type. These programs will also probably be 
used to reduce the OSU Institute of Polar Studies Greenland Doppler data, 
and possibly to provide checks on the geometric Doppler Solutions now 
being made here. 
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